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The dynamical mechanisms that can stabilize the coexistence of species (or strategies) are of substantial in- 
terest for the maintenance of biodiversity and in sociobehavioural dynamics. We investigate the mean extinction 
time in the coevolutionary dynamics of three cyclically invading strategies for different evolutionary processes 
on various classes of complex networks, including random graphs, scale-free and small world networks. We find 
that scale-free and random graphs lead to a strong stabilization of coexistence both for the Moran process and 
the Local Update process. The stabilization is of an order of magnitude stronger compared to a lattice topology, 
and is mainly caused by the degree heterogeneity of the graph. However, evolutionary processes on graphs can 
be defined in many variants, and we show that in a process using effective payoffs the effect of the network 
topology can be completely reversed. Thus, stabilization of coexistence depends on both network geometry and 
underlying evolutionary process. 

PACS numbers: 87.23.-n, 02.50.Ey, 89.65.-s 



I. INTRODUCTION 

Dynamical mechanisms that stabilize cooperation have in- 
trigued scientists from different fields for many decades, 
and evolutionary game theory has developed from a coining 
metaphor to a well-established research area with a wide range 
of applications from biology to behaviour UjtSJ. An analo- 
gous question addresses the dynamical mechanisms that can 
stabilize coexistence of (biological) species or (behavioral) 
strategies [6J. In a biological context: how is biodiversity 
sustained? Here, cyclic coevolution has been suggested as 
to support the coexistence of strategies f?'-^. But it already 
has been reported that cyclic dominance alone is not enough 
to stabilize coexistence of strategies ||8l[l0t]. While space (or, 
in ecological setting, niches) provides a means of stabilizing 
biodiversity, it is not the only mechanism to be taken into ac- 
count. One additional way of stabilization is the introduction 
of a non-zero sum game, which also can stabilize finite but 
sufficiently large populations ifTlll and results in an exponen- 
tial scaling of the mean first extinction time 

In general, the stability of the coexistence fixed point de- 
pends on several instances: firstly, on the payoff matrix 
13, [Hh, secondly, on the population size and the underlying 
evolutionary process lfT3ill4ll . and thirdly, on the spatial struc- 
ture of the population which can assume interaction topolo- 
gies with the structure of a lattice or more general graphs. A 
spatial discretization of coevolutionary dynamics can lead to 
a stabilization of the coexistence of game theoretical strate- 
gies, compared to well mixed populations I7|,l8il- This cor- 
responds to the effect of spatial discretization in 2-strategy 
games, where - despite exceptions II 1511 - cooperative strate- 
gies can be stabilized lfT6l[T7[| . In this direction, many works 
deal with structured versions of the prisoners dilemma ifTsll . 
where it has been shown that cooperators can coexist with 
defectors even in parameter regions where cooperation will 
never be observed in mixed populations lfl9l - l26ll . In con- 
trast to these two-strategy evolutionary games, here we ana- 
lyze the widely known rock-paper-scissors game (RPS) which 
includes three strategies with cyclic dominance. Such three- 
strategy models including cyclic dominance are well-known 



in game theory and related fields Il8ll9l l27l - l29ll . In a spatially 
extended system, the coexistence of the three cyclically invad- 
ing strategies is stabilized, even in cases where it is not stabi- 
lized in the well-mixed system |7]. - But how strong is this 
stabilization? Are there differences between the stabilization 
effect on different graphs and for different processes? How 
does the population size influence the stabilization? To an- 
swer these questions, in this paper we examine three different 
update mechanisms - the frequency-dependent Moran process 
isollsill . the local update process lll3ll . and a process adapted 
from Szolnoki, Perc, and Danku [22] - and a broad range of 
complex network types as underlying structures, and, for each 
network type, analyze the mean extinction time {MET), i.e. 
the average time until one of the three strategies has gone ex- 
tinct. 



II. THE ROCK-PAPER-SCISSORS GAME 

We consider a population of N individuals or agents. N is 
assumed to be constant all the time (although this common as- 
sumption is an approximation, and additional effects emerge, 
e.g., in growing populations fss']). Each individual is placed 
on one vertex of a graph containing N vertices and can take 
on one of the three strategies 'rock', 'paper', or 'scissors'. As 
played on schoolyards, 'rock' looses against 'paper' but wins 
against 'scissors', and cyclically permutated. In game theory, 
the 'interaction kernel' quantifying the outcome of agent col- 
lisions is cast into a payoff matrix, which reads for RPS 



P = 




(1) 



and we assume that we have a zero-sum RPS game. We now 
can find the total payoff for a certain individual i playing strat- 
egy S{i) by adding all payoffs that i gets from neighbored 
nodes. With L as the adjacency matrix (whose elements Li^j 
are 1 if two vertices i and j are connected by a link, and oth- 
erwise), we can compute the total payoff tt; for an individual 
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i by 
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Assume the strategy S{i) to be a. Let us use the notation /3 
for the strategy dominated by a, and 7 for the strategy that 
dominates a, with {0, 1, 2} as numerical values that represent 
these strategies. Then using the special structure of our payoff 
matrix, this term can be simplified to 



the case of zero-sum RPS we have analyzed here (although 
this still holds for the expectation value (tt)). Note also that 
the original definition of the Moran process is slightly differ- 
ent because one should choose an agent i for reproduction at 
random proportional to (f)M- That means we are ensured that 
there is a real update event in every update step, which is not 
the case in our definition. Our definition is numerical cheaper 
because it is not necessary to compute the payoffs for all of 
the vertices in every time step, and it influences the time scale 
only by a constant factor. 



N 



(3) 



This total payoff can be modified for normalization issues, as 
we will describe in the next section, because the way of nor- 
malization varies in the different update mechanisms, defin- 
ing rules after which strategy changes occur, depending on 
TTj. Note that this total payoff can vary relevantly from the 
case of unstructured populations. For example, consider the 
case of a node i playing strategy a. Let j be the only vertex i 
is connected to, and let j be the only vertex that plays strategy 
7 which may dominate a so that i gets a total payoff tt^ = — 1 
although in the whole population there is only one agent that 
dominates strategy a. This single agent playing strategy 7 
would only have negligible influence on the payoff of i in a 
well mixed population, namely, reduce it by 1/{N — 1). 



B. Local update process on networks 

As the availability of the global information (tt) to all in- 
dividuals may often be unrealistic, one should also consider 
processes where only locally available information is of rel- 
evance to the dynamics. In the extreme case, only randomly 
selected pairs of nodes compete, as in the class of local com- 
parison processes ll24l [32l - r361 . From this class, we have inves- 
tigated the local update process ifTsIl as well. Here we, again, 
normalize the payoff of a vertex i by dividing it by the number 
of its neighbors. Then we again choose two finked vertices i 
and j. With probability 
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j changes its strategy to the strategy of i. Here, /S.-Kmax is 
the greatest possible payoff difference, that is the difference 
between the largest and the smallest entry in the payoff matrix. 



III. UPDATE MECHANISMS 

A. The frequency-dependent Moran process on networks 

For the first update rule we have adapted the well known 
frequency-dependent Isill Moran process llsoll . For normal- 
ization, we divide the total payoff tt^ by k{i) (being the vertex 
degree, number of neighbors of i). So the payoff reads 



N 
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C. SPD-process 

As a third update rule we have modified a process intro- 
duced by Szolnoki, Perc, and Danku |22] which one could 
think of as a simplified version of the local update process at 
the first view. We choose a vertex i per random and another 
vertex j which is connected to i by a link also per random. 
In contrast to the Moran and the local update process, both 
payoffs are modified according to the rule 



(8) 



Now we choose a vertex and another vertex j, for which (j„ ^^^^ a = we again have the payoffs Hke in the other 

two processes.) Now, only if tFJ > 717, j reproduces with 
probability 



I/j , = 1 holds, at random. With probability 



,S(i)^S(») 



1 1 — CJ + UJHi 

2 1 -w + w(7r) 



(5) 



^SPD 



= (ttj - TTi)/ max(7rj - iTi) 



(9) 



the individual on the vertex i reproduces, and its strategy re- 
places the strategy of the player on the vertex j. Here a; is the 
strength of selection and 



N 

(tt) = ^iij7rj/fc(i) 



(6) 



proportional to the difference of the payoffs, and the strategy 
of vertex j replaces the strategy of vertex i. The maximum 
in this term has to be understood as the the maximum of the 
possible payoff difference that two vertices with degrees k{i) 
and k{j) can have. This ensures the probability to always be in 
the range of and 1. One can easily compute this maximum, 
so we can achieve the following: 



is the mean fitness of the relevant subpopulation, in our case, 
all neighbors of i. Note that in contrast to the Moran pro- 
cess in well mixed populations, (tt) needs not to be even for 



Here fc„ 



2(1 + a(A:„„^ - 1)) 
is the maximum of k{i) and k{i). 



(10) 



(a) (b) (c) (d) (e) 

FIG. 1: Realizations of five different types of networks with an average vertex degree of fc = 4 and a total number of vertices of 
N — 25. (a) Square lattice with periodic constraints. The half-links are clear illustrations of the periodic constraints, (b) 
Random regular graph. Each vertex has the same degree k. (c) A connected ER-random graph, (d) Scalefree network after 
Barabasi and Albert. The nodes placed on the square in the middle are the fco starting nodes, (e) Smallworld network after 

Watts and Strogatz with p ~ 0.25. 



IV. INVESTIGATED NETWORKS 

All networks used as an underlying topology for such a 
game should fulfill some common properties. They should 
be connected, so that we do not have some parts of the popu- 
lation not being linked to the rest of the population. Here, we 
always consider undirected and unweighted (binary) graphs. 
To avoid self-interactions of the agents, self-loops should be 
excluded. A schematic illustration of the network types used 
in this paper is shown in Fig.[T] 



The graphs we have investigated are constrained in two ways. 
On the one hand, they must be connected, what means that 
there is at least one possibly path from each node to any other 
so that we don't have two or more groups of not interrelated 
agents. On the other hand, we restrict ourselves to networks 
with a mean vertex degree of 4 or 8 for better comparability 
with the regular lattice. Because we will need only connected 
graphs with a certain mean degree k, and because we would 
find these networks only with some probability even by 
adjusting p with respect to the number of nodes, we have used 
the following algorithm: 



A. Regular lattice 

The first step of complexity of placing vertices spatially - 
beyond the linear chain - is a regular square (or rectangular) 
lattice. Typically, each node is connected with 4 or 8 near- 
est neighbors. To eliminate boundary effects, we use periodic 
boundary conditions. This is a very simple, but in many cases 
not realistic model, so we have analyzed some other specific 
types of models that will be presented in the following sub- 
sections. For a better comparability with the regular square 
lattice with natural connectivity of 4 or 8, we have restricted 
ourselves to graphs with an average vertex degree of 4 and 8. 



B. Erdos-Renyi random graph 

The most frequently used model of real, not lattice-like 
graphs in the past decades was the Erdos-Renyi (ER) random 
graph [37]. Here one starts with a graph with TV nodes and 
connects each pair of nodes with probability p. In this model 
the probability of finding a node with degree k follows a Pois- 
son distribution 



where 



P{k) 



N -1 



\\k 
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(11) 



(12) 



1. Start with a graph G with N vertices and no edges. 

2. Link randomly chosen vertex pairs if they are not al- 
ready linked until the total number of edges equals Nk. 
Now the mean vertex degree will equal k. 

3. If the graph is connected, accept it, otherwise return to 
step (1). 



C. Random regular graph 

A compromise between the random graph and the regular 
lattice is the random regular graph with the same degree 
k for all vertices but a much shorter average shortest path 
length than in a lattice. For generating the random regular 
graphs we have used an algorithm proposed by Steger and 
Wormald llssll . This algorithm ensures to give connected 
random graphs with equal degree k for all vertices: 

1. Start with a graph G with N nodes (1, 2, n) and no 
edges. 

2. Repeat the following until the set S is empty: Let S be 
the set of vertex pairs {w, w} of G that are not connected 
by an edge yet with both having at most degree k — 1. 
Then choose a specific pair {u, v} out of S with prob- 
ability proportional to [k — k{u)) {k — k{v)), where 
k{u) and k{v) are the degrees of u and v in the graph 



4 



generated up to now. Add the edge to G. Delete the pair 
{u, v} from S. 

3. If G is /c-regular (that means connected and of the same 
degree k for all vertices), accept it, otherwise return to 
step (1). 

The last step with the return statement is necessary because 
unfortunately the algorithm does not exclude the possibility of 
receiving disaggregated graphs (which occurs rarely) or that 
the last connection possibihties are only self-loops. 

D. Barabasi-Albert scalefree network 

Many real networks, primarily huge ones, are described by 
the ER-model in an insufficient way. In many of these net- 
works, for instance the internet JUt], the WWW [40^1411], the 
network of power nodes in the west of the United States ll42ll . 
the network of scientific collaborations (a^, or the metabolic 
network of yeast ll44ll one observes a degree distribution that 
(approximately) follows a power law, 

P{k)o^k-^. (13) 

Barabasi and Albert ll45i1 have proposed the following algo- 
rithm of modelling complex real networks: 

1. Start with a small number fco vertices that are all con- 
nected to each other 

2. Add another node and link it to an already present ver- 
tex with probability proportional to the degree of this 
vertex. Add further links in the same fashion until the 
new vertex has a degree fc < fco- 

3. Repeat step (2) (until you have reached some value of 
N you want). 

In contrast to the ER-model for random graphs or the like- 
wise frequently used Watts-Strogatz-Modell (i^, see next 
subsection) this model incorporates two main aspects of real 
networks, growth and preferential attachment. Namely, most 
real networks do not have a fixed size but grow like the WWW 
ll40i |41[1 . The WWW is also an example for the preferential 
attachment: When creating a new website, one will typically 
try to link it to well known and popular sites to gain more at- 
tention for ones own site. Because of that, the Barabasi-Albert 
network has a high density of so-called hubs. These are ver- 
tices that are connected with many other nodes with relatively 
small vertex degree. 

E. Small- World-Network after Watts and Strogatz 

Another frequently used network model is the Small-World 
network. The existence of long range connections with si- 
multaneous appearance of a high local cluster coefficient is 
a phenomenon often obtained in real networks [46].). The 
notion integrates some slightly different types of networks. 



Here we have used the original model of Watts and Strogatz 
ll42ll . In their model one starts with a ckcle of nodes. All of 
them are connected with fc nearest neighbors (fc/2 in the one 
direction and fc/2 in the other). Afterwards, one starts with 
rewiring the links, beginning with an arbitrary vertex i on the 
circle. With a certain probability p, i is disconnected from its 
nearest neighbor in clockwise direction and is connected to a 
randomly chosen vertex j that is not already linked to i. In 
this way one goes on in clockwise direction until all links to 
nearest neighbors have been checked for rewiring. Then one 
goes on with the second nearest neighbor, and so on, until all 
links have been checked for rewiring. 

For p — > one achieves obviously some kind of regular 
lattice (that is structured in a different way that the one we 
have used). In contrast to that, for p 1 one achieves a 
graph equivalent to the ER-random graph. But for interme- 
diate values of p one gets networks with high local cluster 
coefficients like a regular graph (or lattice), but also a small 
average shortest path length like a random graph |47J. For the 
sake of clarity, we have usedp = 0.25 in all cases. 

F. Network with Gaussian degree distribution 

As a sixth type of networks, we have used a graph with a 
Gaussian degree distribution. This type is a compromise be- 
tween the ER-random graph and the random regular graph. To 
generate these networks we have modified the algorithm used 
for the random regular graph by making the designated vertex 
degree fc dependent of i, i — 1, 2, n. For these fc(i) we 
have used a Gaussian distribution with a standard deviation of 
fc/3, where fc is the designated average vertex degree. By this 
choice of the standard deviation it is ensured that almost all 
designated degrees k{i) are positive. We have ignored the few 
occurring non-positive (negative or 0) degrees by choosing a 
new designated degree in the case of such a choice. 

G. Uncorrelated scalefree graph 

To analyze possible effects of the vertex degree distribution 
in scalefree graphs after Barabasi and Albert on the dynam- 
ics of our games, we have investigated the dynamics on un- 
correlated scalefree graphs as well. For the construction of 
such networks we have used the algorithm described in the 
preceding subsection (IIV Fb with a power law degree distri- 
bution identical with the distribution of vertex degrees of the 
Barabasi-Albert scalefree graph instead of a Gaussian distri- 
bution. To avoid confusion, we will refer to this graph short 
as uncorrelated. If we refer to scalefree graphs (without men- 
tioning anything about the correlation) we always consider the 
Barabasi-Albert network. 



V. MEAN EXTINCTION TIMES 

We have carried out extensive simulations to obtain the 
mean time until one of the three strategies has gone extinct 
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when starting with each of the strategies having the same fre- 
quency (this requires that N is an integer muhiple of 
three). After the first strategy has vanished, the game lacks 
its cyclic dominance, and we return to the two strategy case 
which is well known and will end up in a much shorter time in 
a single-strategy steady-state because one of the two remain- 
ing strategies definitely dominates the second. Note that we 
measure time in units of single update steps instead of Monte 
Carlo steps, where on average each agent has been updated in 
one time step, so there is a factor of N between these two time 
scales. 



A. Mean extinction time in well mixed populations 

First, let us consider the case of well mixed populations for 
comparison, which is identical to the case of fully connected 
graphs (that means each vertex is directly linked to any of the 

— 1 others). We have to distinguish between two cases: the 
Moran and the local update process on the one hand, and the 
SPD-process on the other hand. 



1. Moran and local update process in mixed populations 

In well-mixed populations, the Moran and local update pro- 
cess have been shown in |13] to have the adjusted replicator 
equation [48] and the standard replicator equation [3] as limits 
as N goes to infinity. They both have a neutrally stable (for the 
zero-sum RPS) fixed point in the state space of strategy densi- 
ties at -J, -j). For any other point inside of the simplex, 
both replicator equations predict neutrally stable oscillations 
around this fixed point, but in the case of finite populations 
the population can run out of the fixed point and the orbits 
around it because of stochastical fluctuations, and end up on 
the edge of the simplex S3 which is the boundary of the state 
space. A point on the edge of this simplex corresponds to the 
extinction of at least one of the three strategies (there are three 
points on the boundary corresponding to the survival of only 
one strategy). So we expect the MET of this system to tend to 
infinity as A^ cx). But at this point there are still two open 
questions: What is the exact scaling dependence of the MET 
on A^? And how does the selection strength uj influence this 
dependence? 

Our simulations give clear and short answers for both ques- 
tions: In the zero-sum RPS game, the MET is proportional 
to A^^ (in single birth-death update steps, corresponding to A^ 
Monte Carlo steps), and independent of uj (both is different in 
the non zero-sum RPS case). For both processes, Moran and 
local update process as well, we have found the mean extinc- 
tion time text to scale as lfT2ll 

(ie.t) = 0.54(±0.02)7V2. (14) 



2. SPD-process in mixed populations 

As this process has never been defined for unstructured 
populations let us use the complete graph as a model for the 
mixed population for simplicity, so that we do not have the 
risk of making any alterations to this update rule if re-defining 
it. Now we find that, when starting with -j, -j), each 
node has as many neighbors of the dominating strategy as of 
the dominated one, namely A^/3. Because of the structure of 
the payoff matrix the fitness is identical to zero for every strat- 
egy and node. So the condition that the second chosen vertex 
has a greater fitness than the first is never fulfilled, and one 
will never observe a change in the strategy densities. For this 
reason the MET becomes infinite even in finite populations (at 
least when starting in the point ( ? ; ) ; even when start- 
ing in some neighboured point in the state space we would 
only achieve some conditional mean extinction time). 



B. Moran process on networks 

From the results of our simulations (see Fig. |2|, there are 
two immediate main observations: First, for some networks 
biodiversity is greatly stabilized as the MET increases expo- 
nentially with A^ , while for other networks biodiversity is only 
slightly stabilized or even destabilized. Second, the MET is 
no longer independent of the selection strength w even in the 
zero-sum RPS case that we have considered here. Stabilizing 
effects of the underlying structure are more distinct for greater 

CJ. 

Great stabilization of the biodiversity is only observed for 
networks with a heterogenous degree distribution and strong 
selection, while for networks with homogenous degree distri- 
bution we have no such strong stabilization. In fact, the net- 
works with homogenous degree distribution can even lead to 
a slight destabilization of biodiversity for strong selection, as 
the MET grows slower than A^^, as we can see for the random 
regular graph. Even in the cases of the small world and the 
Gaussian network, which both have a small heterogeneity in 
their degree distribution, we have such a destabilization. Only 
for the square lattice we have a stabilization (although the lat- 
tice has a homogenous degree distribution, of course), but the 
stabilization is much weaker than in case of the random graph 
or the two scalefree networks. For the stabilization on the lat- 
tice, the large average shortest path length might play a role 
(which is prop ortional to \/N instead of log A^ as for random 
networks ll49ll V 

For small w, there is not a great change in the behaviour of 
the MET on all networks compared to unstructured popula- 
tions (Fig.|2]i. 

Also, the correlation strength of the network degrees has 
only minor effect of the stabilization of biodiversity. The MET 
for an uncorrected scalefree graph is, in all cases considered 
here, somewhat smaller than for the Barabasi-Albert scale- 
free graph, but in comparison to the other effects, mainly the 
enormous stabiUzation impact of scalefree networks, this is 
negligible. 
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Lattice Random 




200 400 600 SOO 1000 1200 

FIG. 2: (Color online.) Mean extinction times of the Moran process on networks. The plots are semi-logarithmic over N, and 
the MET is devided by N'^ for a better comparison with the result of well mixed populations and to better discern stabilizing 
effects. Red circles: k — A, uj ~ 0.05, orange squares: k = A,uj = 0.40, green rhombi: fc = 8, w = 0.05, blue triangles: k — 8, 
u! — 0.40, black points: well mixed population. Red circles on the upper edge of a plot in connection with a curve of a different 
colour mark are only marks that the curve is going on further, but it isn't depicted for a better presentation. Averages are taken 

over a sample of 100 runs (N > 400) or 1000 runs (N = 1000). 
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C. Local update process on networks 

One of our main points of interest was whether the evolu- 
tionary process itself influences the MET in the case of a net- 
worked population. The corresponding results for the Local 
Update are shown in Fig. [3] There is no qualitative difference 
between the results in the Moran process and in the Local 
Update process. Same as for the Moran process, heteroge- 
nous degree distributions cause an enormous stabilization of 
the biodiversity for strong selection (the MET grows exponen- 
tially with N), while homogenous degree distributions do not. 
The MET is now dependent on the selection strength as well, 
and for small lu it is proportional to N'^ as in the neutral case. 
The only difference is that stabilizing effects are not as strong 
as in the Moran process for the same or even a bit larger uj. 



D. SPD-process on networks 

The SPD-process behaves completely different, as shown 
in Fig. m We recall that for this update rule the quite gen- 
eral observation that spatial discretization can lead to a stabi- 
lization of biodiversity does not hold. Even ignoring this and 
taking a MET scaling for large N like the N'^ scaling from 
Moran process and Local Update as reference, for all stud- 
ied graphs except the lattice we have a destabilization, and 
the MET grows slower than iV^ with only small differences 
between the different networks. Only in the lattice case the 
MET grows exponentially with N. 

Changing the parameter a has minor influence on the prop- 
erties of the MET than the selection strength u in the Moran 
and the Local Update process. For degree homogenous net- 
works, the MET is even independent of a. For networks with a 
strong heterogeneity in the degree distribution the MET grows 
with a for the same value of TV. 

There is also some influence of the correlation of network 
degrees. For the uncorrected scalefree network, the influ- 
ence of the parameter a on the MET is not as great as for 
the Barabasi-Albert network, and the METs for the same N, 
but different a and h are closer together But compared to the 
other effects, this is again negligible. 



VI. THEORETICAL APPROACH 
A. Moran process 

1. General case 

For a version of the RPS game with constant reaction 
rates it has already been shown that networks with heteroge- 
nous degree distribution effect a stabilization of biodiversity 
ifsoll . This corresponds with our simulation results so we have 
adapted the approach from [5 0.1 to underpin the shown be- 
haviour theoretically. 

First let us define the probability density of the strategy a 
on all nodes with degree k as Pa,k: wherea = (0, 1, 2). Here 



we again use the notation a for an arbitrary strategy, (3 for the 
strategy dominated by a, and 7 for the dominating strategy. 
Then the probability da' that a neighbored node of an arbitrary 
vertex plays strategy a' reads 



9a' = ^ k'pk'Pa'/{k). 



(15) 



Here is the probability of finding a vertex with degree k. 
Dividing by the average vertex degree, (k) = J^k ^Pk, as- 
sures the standardization. Further, by the conservation of the 
total density one of the quantities pa.k and of 9a is given so 
that we can eliminate p2.k and 6*2 using po.fe and pi^k as vari- 
ables, or 9q and 9i, respectively 



P2,k 
O2 



PO,k — Pl,k 
9o — 9i 



(16) 
(17) 



The average payoff of a player of strategy a on a node with 
degree k is then given by 



(18) 



We find that the average payoff is independent of the vertex 
degree k so we can drop the index k, ■Ka,k = tTq,. With this we 
can compute the average payoff of neighbours of a vertex 

, , J2l=oJ2k'^o,^ji^Pa.k'k 
= k 

a=0 k' ^ ' 

2 

Q=0 k' 
2 

- Y,{9p-9^)9a. (19) 



By using k'pk' Pa' / {k) instead of pk' we consider that is more 
likely to have a vertex with a high degree as a neighbour than 
one with a smaller degree if both are equally frequent. Thus 
we find the probability for a player with strategy a sitting on 
a node with degree k to reproduce as 



1 1 — U! + CJTTq, 
Wa,k = T;PkPa,k- ; ^ 

2 1 — cj + uj{n) 



(20) 



and the rate with that a player of strategy a sitting on a vertex 
with degree k conveys his strategy to a node with degree k', 
which has played strategy a' before, reads 



Pk'Pa',k' 
(fc) 

1 Pk'Pa',k 

2 (fc) 



Wa,k 



-PkPa,k 



(21) 



1 — UJ + UJTTa 

1 — w + w(7r) 



For Pk — 5k,ko (Sk,ko is the Kronecker symbol) this reduces 
to 



Pa Pa 



1 — UJ + UJTTa 

I — uj + uj{tt) 



(22) 



8 




FIG. 3: (Color online.) Mean extinction times of the local update process on networks. The plots are semi-logarithmic over N, 
and the MET is devided by for a better comparison with the result of well mixed populations and to better discern 
stabilizing effects. Red circles: k = A, lu = 0.05, orange squares: k — 4, uj — 0.50, green rhombi: /c = 8, cj = 0.05, blue 
triangles: k — 8, uj = 0.50, black points: well mixed population. Red circles on the upper edge of a plot in connection with a 
curve of a different colour marks that the curve is going on further, but it isn't depicted for a better presentation. 
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FIG. 4: (Color online.) Mean extinction times of the SPD-process on networks. The plots are semi-logarithmic over N, and the 
MET is devided by iV^ for a better comparison with the result of the other processes. Red circles: k = 4, a — 0.0, orange 
squares: k — 4, a ^ 1.0, green rhombi: k — 8, a = 0.0, blue triangles: k = 8, a — 1.0. 
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Here we have omitted any fc-indices because there is only 
one degree left. The payoffs reduce to T^a — Pp — P-y and 
("■) = X]q=o '^aPa^ ™d one obeys the hopping rates for the 
Moran process in well mixed populations Ml 111 independent of 
ko. This rate describes how likely it is that a player with strat- 
egy a' is replaced by a copy of another player with strategy a 
in a single time step. 

Let us now consider a heterogenous degree distribution. In 
the simplest case, such a distribution consists of only two de- 
grees ki and k2 with frequencies Pi ~: p and P2 = ^ — Pi = 
1 — p, respectively. The replicator equation that we obtain 
from the master equation in the limit N oo for the chang- 
ing of the frequency of the strategy on nodes with degree ki 
reads (for the methodology see e.g. lfl3[l5lll52ll ): 



Pom 



->-o 



ki .ki 
^k2,ki 



-'-ki,ki 



rTiO->2 



(23) 



Inserting our special degree distribution, after a longish but 
straightforward calculation we find 



Pom 



Pom 



kip^ 



2(fc)(r+(7r)) 
kip{l-p) 



(24) 



(r (po,fc2 ~ Pom) 



2(fc)(r + (^)) 

+PoM'^o - Po,fci(7r)2) . 

Here F — is the background fitness. (7r)i = -koPom + 
T^ipiM + ^2P2,ki is the mean fitness of all players on vertices 
with degree ki, and (7r)2 correspondingly for fc2. Then one 
can compute the other equations in an analogical way or by 
cyclically permuting the indices. For each vertex degree, one 
of the equations of motion needs not to be taken into account 
by the constraint po.fc; + Pi.k^ + P2.ki = 1- For ~ ^2 or 
alternatively p — I this flattens to the well known adjusted 
replicator equation in the case of well mixed populations, as 
expected iflslfTill . Likewise, the fixed point of this replicator 
equation at (i, i, i, i) does not alter, we now just have two 
densities for every strategy. But now the velocity of the reac- 
tion is directly dependent of oj because it is no longer possible 
to absorbe the strength of selection by a dynamical rescaling 
of time, as it is possible for well mixed populations. 

But what about the stability of this inner fixed point? 
Within a linear stability analysis, let us linearize the replicator 
equation around its fixed point 



p = Ap, 



(25) 



withp= (poM'PiM^PoM^PiM) and ^ being the Jacobian 
of the system with ^i, — i i i\. In our case, we 

J apj 'P— {3,3,-3,3 ) 

find A as 



/ fci(p-i)p 

' 2<fe> 
, 2 2 

fc2(p-l)P 
2(fc) 

fclfc2(p-l)p 



_2_iP_ 

6(fcpr 
fci(p-i)p 

2{k) 
fclfc2(p-l)P 

fc2(p-l)P 



fcl(p-l)p fclfc2(p-l)p \ 

2{k) 6(fc)^r ' 



fclfc2(p-l)p 

6(fc>^r 
fc2(p-i)p 

fci(p-i)" 

6(fc>^r 



fci(p~i)p 

2(fc) 

fci(p-i)" 

6(fc)^r 
fc2(p-i)p 

2(fc> 



To reason about the stability we need to know the signs of the 
real parts of the eigenvalues A,; of this matrix, but for these 
eigenvalues we derive quite complicated terms, the signs of 
which cannot be easily designated (but see AppendixlAli. The 
Routh-Hurwitz criterion as an alternative solution does not 
simplyify the situation considerably. 

Approximating an upper bound of the real parts of all eigen- 
values, it is possible to compute a region in the parameter 
space where the inner fixed point is guaranteed to be asymp- 
totically stable, namely, for 



0.5 <p < 



(27) 



For a given value of p this determines the ratio of kmax and 
kmin, or, the other way round, if we want to build up a net- 
work with chosen values for k^ax and fc„ii„, this inequality 
tells us what frequencies of both vertex degrees are necessary 
to stabilize biodiversity. So, for example, if we choose ki — 4 
and both vertex degrees should be equally frequent, k2 must 
be at least 7 to fulfill this inequality. If a higher frequency of 
ki is desired, k2 must even be larger. 

But as this is just a rough approximation, by inserting nu- 
merical values we have found that almost all eigenvalues are 
< 0, except for the trivial cases fci = k2, p — or p ~ 1. So 
for all values except these mentioned cases, we have found the 
fixed point to be asymptotically stable. Here, we have found 
numerically that the real parts of the eigenvalues become more 
negative as the difference of ki and k2 grows in absolute size, 
and both vertex degrees become equally likely to be found in 
the network. This explains as well why the stabilization of 
biodiversity is greater in the case of the scalefree graphs than 
for the ER-random graph, or it is in both cases greater than 
for the network with a gaussian degree distribution, which has 
only as small heterogeneity in its degree distribution. 



2. Small selection limit 

As we have seen in the previous section, for small lu the 
MET is almost proportional to N"^ as in the well mixed case, 
and depends only marginally on the underlying network struc- 
ture. We can intuitively understand this as follows: If u) is 
small, than every transition probability is « i, and the corre- 



sponding replicator equation in our system reads 

Pa,k = 



(28) 



(26) 



for all values of a and k. Hence, each point of the state space 
is a neutrally stable fixed point in the limit — ^ oo. Therefore 
there should be no stabilization impact compared to the well 
mixed case, as there each point in the simplex belongs to a 
neutrally stable limit cycle, so in both cases a perturbation is 
neither expected to decay nor to grow in time. 



B. Local update process 

Following the same scheme as for the Moran process, one 
can derive similar results (see Appendix IB] for details). As in 
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or 



(33) 



PO k ?P0 k2 I'^sP^'-'^i^^ly- account for the fact that ttj > tt^ must 



FIG. 5: (Color online.) A plot of a numerical solution of the 
replicator equation for the local update process on a network 
with only two different vertex degrees (eq. lB2l i with p = 0.5, 

ki=4,k2^8,uj = 0.4, po,fci (0) = 0.65, pi,k, (0) = 0.3, 
Pom (0) — 0.3, pi^k2 (0) = 0.1, red for the density of strategy 
on nodes with degree ki, black for the corresponding value 
on vertices with degree k2 



the Moran process, a degree heterogeneity leads to a stabiliza- 
tion of biodiversity, but again the eigenvalues of the Jacobian 
are too byzantine for catching the signs of the real parts easily. 



C. SPD-process 

For this update process, we will use a simpler method to 
achieve a theoretical understanding of some of the numerical 
results. In the case of both chosen vertices having the same 
degree, the probabihty for the node i to change its strategy to 
the strategy of j simplifies to 



hold, then the first fraction is < 1 for iTj > and > 1 for 
TTj < 0. Hence, if aside from a all conditions are the same, 
for a = 1 it is more likely that a strategy with a lower payoff 
reproduces, if k{j) < k{i)) than for a = 0. On the other 
hand under these conditions for a ~ 1 it is less likely that a 
strategy with a high payoff reproduces, than for a = 0. In 
the second case, the fraction is < 1 for ttj > and tt^ < 
and > 1 for all other cases. Hence for a = 1 it is more likely 
that a strategy with a high payoff replaces a strategy with a 
smaller payoff from a node with a smaller vertex degree, than 
for a — Q. And for a — 1 under this conditions it is as much 
more likely that a strategy with a smaller payoff reproduces. 
Taken together, for a = 1 strategies with a smaller payoff 
have a bigger chance to reproduce, and they are less likely to 
be replaced by a strategy with a higher payoff, than for a = 0. 
Hence on degree-heterogenous networks the MET should be 
greater for a = 1 than for a = 0, as it is confirmed by the 
simulation results (Fig.|4|. 



(29) 



So for the lattice and the random regular graph, where all ver- 
tex degrees are the same, any dependence of a drops out. Note 
that in this equation and in this whole subsection, when talk- 
ing about TTi we consider the not modified payoffs. For a = 1 
one achieves the same term as in the case of degree homoge- 
nous networks (eq.|29l), independent of the network, while for 
a = we find 



Ws(i)^Sij)\a=0 



Hf) 



or 



W< 



S{i)^S(j)\a=0 



(30) 



(31) 



respectively, depending on which of the vertex degrees is 
greater. If we compare the probabilities for a = 1 and for 
a = by computing the ratio of both, we can analyze the im- 
pact of the parameter a and the underlying network structure 
on the MET. We find 



Ws(t)^SU) \a=l _ k{j){Tri - TTj) 



Ws{i)~^S{j} |a=0 fcOVi - kr, 



(32) 



VII. CONCLUSIONS 



Cyclic coevolutionary dynamics manifests an interesting 
class of dynamical systems that are discussed as models to 
explain long-term stabilization of coexistence of species or 
strategies. In this paper we have shown that, for the coevo- 
lutionary RPS game on complex networks, there are consider- 
able differences in the mean extinction time (MET) between 
different networks, and for different processes as well. The 
simple statement "Spatial discretization stabilizes biodiver- 
sity" must be refined in a more sophisticated way. While for 
example in the Moran and in the local update process degree- 
heterogenous networks preserve biodiversity (the MET grows 
exponentially with the number N of players), in the SPD- 
process they do not, and even degree-homogeneous networks 
may destabilize while the lattice still stabilizes. In summary, 
we have shown that, for the cyclic evolutionary RPS on net- 
works, there is a striking influence on the mean extinction time 
characteristics of the spatial structure and the update rule as 
well. 
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Appendix A: Eigenvalues of the Jacobian in equation ( 



1 < cos(a::) > 1, all real parts are smaller than 



With the help of (e.g.) Mathematica'''^ one can easily com- 
pute the eigenvalues Xi of the Jacobian in equation ( |26] ): 



Ai = - 



1 



i2(fc)2r 

3(A:i + fc2)(fc) {p - l)pr + Vci - 1C2) 
^ {-^{ki{p-lr + klp') 



(Al) 



(A2) 



i2(fc)2r 

+3{ki + k2){k){p - l)pT + Vci - 102) 

i^kiip - I f + ik2{k)pT{p - 1) (A3) 



A4 



+kip{ikip + 3{k){p - l)r) - Vci + ic2) 

(iklip - if + 3k2{k)pr{p ~ 1) 



1 



i2(fc)2r 



(A4) 



+kip{ikip + 3{k){p - l)r) + Vci + ica) 
with the abbreviations 

'2/ in2 , ,2 2 



(A5) 



ci = -{kiip^if + kip^y 

+9{k,+k2f{kf{p-lfp'T' 
C2 = 6{k)ip-l)p{p'kl + k2{p~2)pkl (A6) 
+ki{p'-l) k, + kl{p-lf) r. 



The real parts 5RAi of this eigenvalues are then given by 



nX2 = 



3?A4 = 



fcip2 


/C2P^ kiP k2P 


4(fc) 


^ 4:{k) 4:{k) A{k) 


cos ( 


iarg(ci - ic2)) \/4 + cl 




12(fc)2r 


fcip2 


fep^ kiP k2P 


m 


4{k) 4(fc) 4(fc) 


cos ( 

+ — ^ 


iarg(ci -ica)) -v/c| + cj 




12(fc)2r 


fcip2 




4(fc) 


4(fc) 4(fc) 4(fc) 


cos(i 


arg(ci +ic2)) ^4 + cf 




12(fc)2r 


fcip2 


fep^ fcip k2P 


m 


4{k) A{k) 4(fc) 


cos ( 


iarg(ci +ic2)) i/cj + cl 


i2(fc)2r 



(A7) 



(A8) 



(A9) 



(AlO) 



Here, arg(z) is the argument of the complex number z. The 
sum of the first four terms negative in each case because < 
p < 1 and for this reason p' < p. The problem is the last 
summand, but it is possible to make a rough approximation 



of where all real parts are negative. As 



i2(fc)^r 



> and 



fciP^ fc2P^ _ kiP_ _ k2P_ \/cl + 4 
A{k) 4(fc) 4(fc) 4(fc) 12(fc)2r 



(All) 



Inserting ci and C2 and assuming p > 0.5 without loss of 
generality, one can approximate the fourth root in the last term 
by 

^cl + ci = (36{fc)^(l-p)Vr'(p'fc? + (-2+p) 

2,.3\2 



pklk2 + (-1 + p')fclfc2' + (1 - Pfk'2) 
+ {-{kf{l-pfp'V\k,+k2f 

+ (p^fc? + (i-p)^fci)')')' 

< (36(fc)2(l-p)Vr^(p''fe? + (2 + p) 



pfc?fc2 + (1 +p')fcifci + (1 - pfklf 

+ ((fc)2(i-p)Vr'(fci + fc2)' 

+ [p'ki + ii-pfkify 

< (36fcL.p'r2(fcL.(4p' + 2))^ 

< (36 ■ 144r2/fcl,, + (36r2 + 4)2/ fcl,,) ^ 

= p'fcL. \/36 ■ 144r2 + (36r2 + 4)2 (A12) 

As the have checked numerically, the remaining root divided 
by ^36 • 144r2 + (36r2 + 4)2/r < 9 for cj < 0.45 (which 
is necessary to have mathematically meaningful probabilities 
between and 1), so 



V4Tci ^ 3p2fc: 



2 

'max 



i2(fc)2r 4(fc) ' 

and therefore all real parts are negative if 

kniax 4* kjnin ^ P [ ^?na2: ^ k^nin 4~ 3 



(fc) 



(A13) 



(A14) 



or noted in a simplified way and remembering the assumption 
at the beginning of this approximation. 



0.5 <p < 



(A15) 



Appendix B: Theoretical approach: local update process on 
networks 

Similar as in the Moran process, in the local update process 
on networks the rate with which vertices with degree k and 
strategy a carry this strategy over to nodes with degree k' and 
strategy a' is given by 

ya'^-^a ^ Pk'Pkk'pa,kPa',k' |^ 1 ^ I ^ T^g -TTg' ^ ^^^^ 



(k) 



2 2 Att,, 
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which does not vary from the results in well mixed popu- 
lations in the case of homogenous degree distribution. The 
replicator equation for the model system of a heterogenous 
degree distribution with only two vertex degrees reads 



Pom (tto - (7r)i) 



+ 



(k) 
kip{l-p) 



(B2) 



(po,fe2 (l + V'(7i"o(7r)i)) 



2{k) 

Pom (1 + V ((7I")2 - tto))) 



with tp = u/A-jTmax- The other equations are obtained anal- 
ogously. Here again the position of the fixed point is not 
changed compared to the well mixed case. But again it is no 
longer possible to absorb the selection strength a; by a rescal- 
ing of time. 

A linear stability analysis leads to the matrix A 



I fci(p-i)p 

fc2(p-l)p 
2<fe> 

felfc2(p-l)pV' 



fcl(p-l)p 

fclfc2 (P — 1)PV' 
3(fcp 
fc2(p-l)p 
2(fc) 



fcl(p-l)P 
1{k) 

felfe2(p-l)pV' 

WP- 

fc2(p-l)p 
2(fc) 

fci(p-i)^v. 



fcl(p-l)p 

2<fc> 

3<fcp 
fc2(p-l)p 
2<fc> 



(B3) 



The eigenvalues of this matrix are again given by 
1 



{-3{ki + k2){k){p-l)p 



12(A;)2 

+2i - 1)^ + klp"^) ip + Vci - ic2) 

^2=Y^mh + h2){k)ip-i)p 

~2i {kl{p - If + klp"^) ^ + Vci - ic2) 



{2ikltpp'^ + 3ki{k){p - l)p 



(B4) 



(B5) 



12{ky 



(B6) 



+k2{p - l){3{k)p + 2ik2{p - l)tp) - Vci + ic2) 
Xi =7^77^ (2ifc?Vp^ + 3ki{k){p - l)p (B7) 



12(fc)2 

+k2{p - l)(3(fc)p + 2ik2{p - l)ip) + Vci + ic2) 



with the abbreviations 



2„2 



Ci = 9(fci+fc2)'(fc)"(p-l)V 

-4(fc2(p-l)2+fc?p2)2^2 
C2 = 12(fc)(p-l)(p2fc3+fc2(p-2)pfc2 

The real parts of this eigenvalues are given by 



fcip^ ^ k2p'^ kip k2P 



+■ 



4(fc) 4(fc) 
cos (i arg (ci - ?:c2)) ^cf+cf 



^ k2P' 



12(fc)2 



fclp^ 

'"^ 4(A;) 4(fc) 



fcip k2P 



cos (I arg (ci + ic2 



KA4 = 



A{k) 4(fc) 4(A;) "i^/t/ 

cos (I arg (ci -|- ic2)) -^/c^ +c| 
^ 12(fcp ' 



12(fc)' 
k2P^ kip 



0) y/^FT^ 



k2P 



(B8) 
(B9) 



4(A;) 4(fc> 4(fc) 4(fc) ^ ^ 

cos (i arg (ci - 102)) -|- c| 
12(fc)2 



(BID 



(B12) 



(B13) 



which is very similar to the corresponding result in the Moran 
process (it varies only slightly in ci and C2). Searching for 
an upper bound for the real parts of the eigenvalues, we even 
derive the same inequality as for the Moran process 



0.5 <p< 



(B14) 
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